﻿#include "discrete_edge.hxx"
#include "discrete_face_utils.hxx"
#include <algorithm>
#include <cmath>
#include <iostream>

#define IF_OCCT_ANGLE 1

static int ep_cnt = 0;
static int idx = 0;
double GetSegAngle(double dt, double at, double radius) {
    auto val = 1.0 - dt / radius;
    val = std::min(1.0, std::max(-1.0, val));
    auto angle = 2.0 * acos(val);
    angle = std::min(at * M_PI / 180.0, angle);
    return angle;
}
void similar_triangle_intersection(double x, double prev_x, double adjusted_x, double prev_y, double adjusted_y, double& proportion, double& val) {
    proportion = (x - prev_x) / (adjusted_x - prev_x);
    val = prev_y + (adjusted_y - prev_y) * proportion;
}

/*-----对给定的边对象进行离散处理-----*/
/**
 * EDGE* e：几何曲线的一条边
 * double dt:距离容差
 * double at：角度容差
 * AF_POINT*& out_points：指向AF_POINT对象的指针，存储处理后的几何点链表
 * int& out_cnt：处理后的几何点链表的长度
 */
 // @todo 命名信息可读性
outcome facet_edge(EDGE* e, double dt, double at, AF_POINT*& out_points, int& out_cnt) {
    // 若边指针为空
    if (e == nullptr) return outcome();

    // 如果边的几何信息为空
    if (e->geometry() == nullptr) {
        // 创建一个AF_POINT对象作为头结点
        AF_POINT* head = ACIS_NEW AF_POINT(0, nullptr);
        // 设置参数和位置
        head->set_parameter(e->start_param());
        head->set_position(e->start_pos());
        // 将头结点附加到边上
        head->attach(e);
        out_points = head;
        out_cnt = 1;
        return outcome();
    }

    // 获取边的几何方程
    const curve& dup_curve = e->geometry()->equation();
    int sense = e->sense();
    // 获取边修正后的起始参数和终点参数
    double st = (e->sense()) ? (-e->start_param()) : (e->start_param());
    double ed = (e->sense()) ? (-e->end_param()) : (e->end_param());

    double len = dup_curve.length(st, ed, TRUE);

    // 三维空间中该曲线的长度为0 || 该曲线长度不为0，但是两点过于靠近
    if (is_zero(len) || same_point(st, ed, dup_curve) && fabs(st - ed) < SPAresfit) {
        // @warning 若是有两点在三维空间靠近，但是二维空间远离，需要作为VERTEX而存在，
        // 否则无法通过Brep体有效性的检查,
        AF_POINT* head = ACIS_NEW AF_POINT(0, nullptr);
        head->set_parameter(st);
        head->set_position(e->start_pos());
        head->attach(e);
        out_points = head;
        out_cnt = 1;
        return outcome();
    }

    // 首尾点参数值不接近的情况
    // 是否需要检查容差（根据公式进行分段的椭圆线（包括圆和椭圆）以及螺旋线无需进行迭代容差检查）
    int flag = 0;

    // params.push_back(e->start_param());
    // 新建首尾AF_POINT节点，设置参数和位置，并将边指针附加到头结点上
    AF_POINT* head = ACIS_NEW AF_POINT(0, nullptr);
    // 设置头结点参数和位置
    head->set_parameter(st);
    head->set_position(dup_curve.eval_position(st));
    // 将头结点附加到边上
    head->attach(e);

    AF_POINT* tail = ACIS_NEW AF_POINT();
    // 设置尾节点参数和位置
    tail->set_parameter(ed);
    tail->set_position(dup_curve.eval_position(ed));
    tail->set_prev(head);
    tail->set_next(head);

    head->set_next(tail);
    head->set_prev(tail);

    // 根据度数初始分段
    // @todo 1魔数
    if (dup_curve.type() == straight_type) {  // straight
        // 若是直线，则检查头结点和下一结点的位置是否相同，相同则表示链表中只有一个节点
        if (head->get_position() == head->next(0)->get_position()) {
            head->set_next(head);
            head->set_prev(head);
        }
        // 不需要检查容差
        flag = 1;
    }
    else if (dup_curve.type() == ellipse_type) {
        // ellipse
        // flag = 1;
        // 将曲线转换为椭圆类
        const ellipse* ell = (const ellipse*)&dup_curve;
        // 获取椭圆半径比例
        double radius_ratio = ell->radius_ratio;
        // 半径比例为1说明是圆
        if (radius_ratio == 1.0) {
            // 不需要检查容差
            flag = 1;
            // 调用接口离散圆形曲线
            facet_edge_circle(ell, st, ed, dt, at, sense, head);
        }
        else {
            // 不需要检查容差
            flag = 1;
            // 调用接口离散椭圆曲线
            facet_edge_ellipse(ell, st, ed, dt, at, sense, head);
        }
    }
    else if (dup_curve.type() == helix_type) {  // helix
        // 不需要检查容差
        flag = 1;
        // 将曲线转换为螺旋线类
        const helix* heldef = (const helix*)&dup_curve;
        // 调用接口离散螺旋曲线
        facet_edge_helix(heldef, st, ed, dt, at, sense, head);
    }
    else if (dup_curve.type() == intcurve_type) {  // int-curve
        // 需要检查容差
        // flag = 1;
        // 将曲线转换为插值曲线类
        const intcurve* intc = (const intcurve*)&dup_curve;
        // 调用接口离散插值曲线
        facet_edge_intcurve(intc, st, ed, dt, at, sense, head);
        // if(intdef != nullptr) ACIS_DELETE intdef;
    }

    // flag=1时检查容差
    if (!flag) {
        // 根据距离容差检查离散点，若不满足容差要求，对不满足分段进行二分细分
        CheckDistanceTolerance(e, dt, at, head, head);
        // 根据角度容差检查离散点，若不满足容差要求，对不满足分段进行二分细分
        CheckAngularTolerance(e, at, head);
    }

    CheckSingularities(e, head, dt, at);

    AF_POINT* cur = head;
    AF_POINT* nxt = nullptr;
    AF_POINT* prv = nullptr;
    AF_POINT* dup = nullptr;
    // 循环遍历链表，检查并删除重复的节点
    do {
        // 如果当前节点的位置与下一个节点位置相同
        if (cur->get_position() == cur->next(0)->get_position()) {
            // 如果下一个节点是尾节点
            if (cur->next(0) == tail) {
                // 将下一个节点的下一个节点连接到当前节点，再删除下一个节点
                cur->next(1)->set_next(cur->next(0));
                cur->next(0)->set_prev(cur->next(1));
                dup = cur;
                cur = cur->next(1);
                if (dup != nullptr) ACIS_DELETE dup;
            }
            else {  // 如果下一个节点不是尾节点
                // 将下一个节点的下一个节点和下下个节点分别赋值给nxt和prv，然后将当前节点的下一个节点设为 nxt，删除重复节点
                dup = cur->next(0);
                nxt = dup->next(0);
                prv = dup->next(1);
                cur->set_next(nxt);
                cur->next(0)->set_prev(prv);
                ACIS_DELETE dup;
            }
        }
        else  // 如果当前节点的位置与下一个节点位置不相同
            cur = cur->next(0);
    } while (cur != tail);  // 直到遍历到尾节点结束循环

    cur = head;
    // 检查边界参数域上的点是否需要根据曲面的周期性进行调整
    NewCheckParamBoundary(e, cur);

    // 统计边离散点数量
    int cnt = 0;
    cur = head;
    cnt++;
    cur = cur->next(0);
    while (cur != head) {
        cnt++;
        cur = cur->next(0);
    }
    out_cnt = cnt;
    out_points = head;

    return outcome();
}

/*-----检查边界参数域上的点是否需要根据曲面的周期性进行调整-----*/
/**
 * EDGE* e：几何曲线的一条边
 * AF_POINT* head：指向AF_POINT对象的指针，代表几何点的链表的起始点
 */

 /*-----检查边界参数域上的点是否需要根据曲面的周期性进行调整-----*/
 /**
  * EDGE* e：几何曲线的一条边
  * AF_POINT* head：指向AF_POINT对象的指针，代表几何点的链表的起始点
  */

void NewCheckParamBoundary(EDGE* e, AF_POINT* head) {
    // 获取与边相邻的所有面
    ENTITY_LIST face_list;
    // 将与边e相邻的面添加到face_list中
    api_get_faces(e, face_list);
    auto sense = e->sense();
    auto curv = &e->geometry()->equation();
    double st = (e->sense()) ? (-e->start_param()) : (e->start_param());
    double ed = (e->sense()) ? (-e->end_param()) : (e->end_param());
    auto findBoundaryPoints = [&](const curve* sub_curv, const curve* curv_boundary, double prev_par, double nxt_par) -> curve_curve_int* {
        SPAbox const& box = *(SPAbox*)NULL_REF;
        if (sub_curv == nullptr || curv_boundary == nullptr) return nullptr;
        curve_curve_int* int_res{ nullptr };
        double base_fit = SPAresabs;
        for (int i = 0; i < 3; i++) {
            int_res = int_cur_cur(*sub_curv, *curv_boundary, box, base_fit);
            if (int_res != nullptr) {
                double par = int_res->param1;
                if (curv->periodic()) {
                    par = resetToInterval(par, st, ed, curv->param_period());
                }
                if (is_less_than(par, std::max(prev_par, nxt_par)) && is_greater_than(par, std::min(prev_par, nxt_par))) {
                    int_res->param1 = par;
                    return int_res;
                }
            }
            base_fit *= 10;
        }
        return nullptr;
        };
    auto SolveInterSect = [&](BoundaryAnalyzer& ba, auto GetCord, auto NewPos, double low, double high, bool is_u) {
        AF_POINT* cur = head;
        double period = is_u ? ba.info_.surf.param_period_u() : ba.info_.surf.param_period_v();
        double threshold = period / 2;

        auto curv_low = is_u ? ba.info_.surf.v_param_line(low) : ba.info_.surf.u_param_line(low);
        auto curv_high = is_u ? ba.info_.surf.v_param_line(high) : ba.info_.surf.u_param_line(high);
        while (true) {
            AF_POINT* nxt = cur->next(0);
            // 检查下一个点是否为链表的头部，如果是，则遍历完成，跳出循环
            if (nxt == head) break;
            // 获取前后两点在面上的参数域坐标（二维）
            double prev_par = cur->get_parameter(), nxt_par = nxt->get_parameter();
            SPAposition prev_pos = cur->get_position(), nxt_pos = nxt->get_position();
            SPApar_pos prev_par_pos = ba.info_.surf.param(prev_pos), nxt_par_pos = ba.info_.surf.param(nxt_pos);
            // 可能产生跨界行为
            if (abs(GetCord(prev_par_pos) - GetCord(nxt_par_pos)) > threshold) {
                // 根据前后两点v差值的符号，判断后点需要调整的周期的符号（即，从上界跨越或从下界跨越）
                double diff = GetCord(nxt_par_pos) - GetCord(prev_par_pos);
                double period_diff = period;
                if (diff > 0) period_diff *= -1.0;
                // 计算下一个点在u方向上的调整后的参数坐标
                SPApar_pos adjusted = NewPos(nxt_par_pos, period_diff + GetCord(nxt_par_pos));
                if (ba.SamePointJudgeBy3D(prev_par_pos, adjusted, SPAresabs)) {
                    cur = cur->next(0);
                    continue;
                }
                if (diff > 0 && ((GetCord(prev_par_pos) - low) * (GetCord(adjusted) - low) <= 0)) {  // 与下界相交
                    // 要求点不在界上
                    if (!ba.SamePointJudgeBy3D(prev_par_pos, NewPos(prev_par_pos, low), SPAresabs)) {
                        AF_POINT* down = ACIS_NEW AF_POINT();
                        auto sub_curv = curv->subset({ prev_par, nxt_par });
                        auto int_res = findBoundaryPoints(sub_curv, curv_low, prev_par, nxt_par);
                        ACIS_DELETE sub_curv;
                        if (int_res != nullptr) {
                            AF_POINT* down = ACIS_NEW AF_POINT();
                            down->set_parameter(int_res->param1);
                            down->set_position(int_res->int_point);
                            down->insert(cur, sense);
                        }
                    }
                }
                else if (diff <= 0 && ((GetCord(prev_par_pos) - high) * (GetCord(adjusted) - high) <= 0)) {
                    if (!ba.SamePointJudgeBy3D(prev_par_pos, NewPos(prev_par_pos, high), SPAresabs)) {
                        AF_POINT* down = ACIS_NEW AF_POINT();
                        auto sub_curv = curv->subset({ prev_par, nxt_par });
                        auto int_res = findBoundaryPoints(sub_curv, curv_high, prev_par, nxt_par);
                        ACIS_DELETE sub_curv;
                        if (int_res != nullptr) {
                            AF_POINT* up = ACIS_NEW AF_POINT();
                            up->set_parameter(int_res->param1);
                            up->set_position(int_res->int_point);
                            up->insert(cur, sense);
                        }
                    }
                }
            }
            cur = cur->next(0);
        }
        ACIS_DELETE curv_high;
        ACIS_DELETE curv_low;
        };

    // 开始遍历face_list中的每一个面
    for (FACE* f = (FACE*)face_list.first(); f != NULL; f = (FACE*)face_list.next()) {
        // 如果面的u向和v向均非周期性，跳过此面
        if (!f->geometry()->equation().periodic_u() && !f->geometry()->equation().periodic_v()) continue;  // 面u向与v向上均非周期，跳过

        // 获取面的参数域包围盒、周期以及跨参数边界的阈值
        SPApar_box face_par_box = f->geometry()->equation().param_range();
        double low_u = face_par_box.low().u, low_v = face_par_box.low().v;
        double high_u = face_par_box.high().u, high_v = face_par_box.high().v;
        BoundaryAnalyzer ba(f);
        if (ba.info_.surf.periodic_u()) {  // 面在u向上呈周期性
            // 获取曲面u向上的周期与跨越阈值
            SolveInterSect(
                ba, [](SPApar_pos pos) { return pos.u; }, [](SPApar_pos pos, double cord) { return SPApar_pos(cord, pos.v); }, low_u, high_u, true);
        }
        if (ba.info_.surf.periodic_v()) {  // 面在v向上呈周期性
            // 获取曲面v向上的周期与跨越阈值
            SolveInterSect(
                ba, [](SPApar_pos pos) { return pos.v; }, [](SPApar_pos pos, double cord) { return SPApar_pos(pos.u, cord); }, low_v, high_v, false);
        }
    }
    return;
}

/*-----根据距离容差检查离散点，若不满足容差要求，对不满足分段进行二分细分-----*/
/**
 * EDGE* e：几何曲线的一条边
 * double& dt：距离容差
 * AF_POINT* head：指向AF_POINT对象的指针，代表几何点的链表的起始点
 * AF_POINT* head：指向AF_POINT对象的指针，代表检查区间结束点
 */
void CheckDistanceTolerance(EDGE* e, double dt, double at, AF_POINT* head, AF_POINT* end_pt) {
    // cur指针指向head，用于遍历几何点链表
    AF_POINT* cur = head;
    // 获取边e的几何形状的数学方程
    auto& eq = e->geometry()->equation();

    // 获取边e的生成点的方向sense
    auto sense = e->sense();

    while (cur != nullptr) {
        // 获取当前点cur的下一个点nxt
        AF_POINT* nxt = cur->next(0);
        if (nxt == end_pt) break;

        // 获取端点参数值
        double st = cur->get_parameter(), ed = nxt->get_parameter();

        // 若两端点参数差值较小，不做容差检查
        if (is_equal(st, ed)) {
            cur = cur->next(0);
            continue;
        }
        // 获取起止两点的切线方向 (向量)
        SPAunit_vector vec_p0Tangent = eq.eval_direction(st, TRUE, TRUE);
        SPAunit_vector vec_p1Tangent = eq.eval_direction(ed, TRUE, TRUE);

        auto dir_3d = normalise(nxt->get_position() - cur->get_position());
        if (e->sense()) dir_3d = -dir_3d;
        auto a_1 = acos(std::clamp(vec_p0Tangent % dir_3d, -1.0, 1.0)) / M_PI * 180;
        auto a_2 = acos(std::clamp(vec_p1Tangent % dir_3d, -1.0, 1.0)) / M_PI * 180;
        double w = (a_2 + at / 2) / (a_1 + a_2 + at);
        // 获取两端点的参数域中点
        double mid_t = st + (ed - st) * w;

        // 获取两端点的三维坐标和参数域中点的三维坐标
        SPAposition p0 = cur->get_position();
        SPAposition p1 = nxt->get_position();
        SPAposition mid_p = eq.eval_position(mid_t);

        // 获取两端点连线的中点
        SPAposition mid_3d = p0 + (p1 - p0) * w;

        // 判断连线中点与参数中点映射到三维空间后两点之间的距离
        auto vec = mid_3d - mid_p;

        if (vec.len_sq() > dt * dt) {  // 不满足距离容差，二分区间，迭代检查
            // 创建一个新的AF_POINT类型几何点
            AF_POINT* nw = ACIS_NEW AF_POINT();
            // 设置该点的参数为当前起止段中点参数值
            nw->set_parameter(mid_t);
            // 设置新点的位置
            nw->set_position(mid_p);
            // 将新点插入当前点cur后
            nw->insert(cur, sense);
        }
        else {
            cur = cur->next(0);
        }
    }
}

/*-----根据距离容差检查离散点，若不满足容差要求，对不满足分段进行二分细分-----*/
/**
 * EDGE* e：几何曲线的一条边
 * double& at：角度容差
 * AF_POINT* head：指向AF_POINT对象的指针，代表几何点的链表的起始点
 */
void CheckAngularTolerance(EDGE* e, double& at, AF_POINT* head) {
    // cur指针指向head，用于遍历几何点链表
    AF_POINT* cur = head;
    // 获取边e的几何形状的数学方程
    const curve& eq = e->geometry()->equation();

    REVBIT sense = e->sense();
    // 循环遍历整个几何点链表
    while (cur != nullptr) {
        // 获取当前点cur的下一个点nxt
        AF_POINT* nxt = cur->next(0);
        if (nxt == head) break;

        // 获取当前点和下一个点所对应的参数值
        double st = cur->get_parameter(), ed = nxt->get_parameter();
        // 若两端点参数差值较小，不做容差检查
        if (same_point(st, ed, eq)) {
            cur = cur->next(0);
            continue;
        }

        // 获取起止两点的切线方向 (向量)
        SPAunit_vector vec_p0Tangent = eq.eval_direction(st, TRUE, TRUE);
        SPAunit_vector vec_p1Tangent = eq.eval_direction(ed, TRUE, TRUE);

        auto dir_3d = normalise(nxt->get_position() - cur->get_position());
        if (e->sense()) dir_3d = -dir_3d;
        auto a_1 = acos(std::clamp(vec_p0Tangent % dir_3d, -1.0, 1.0)) / M_PI * 180;
        auto a_2 = acos(std::clamp(vec_p1Tangent % dir_3d, -1.0, 1.0)) / M_PI * 180;
        double w_1 = (a_1 + at) / (a_1 + a_2 + at * 2);
        double w_2 = (a_2 + at) / (a_1 + a_2 + at * 2);
        double mid_t = st * w_1 + ed * w_2;
        // 计算起止两点切线的夹角
        double dot_0_1 = std::clamp(vec_p0Tangent % vec_p1Tangent, -1.0, 1.0);
        double angle_0_1 = acos(dot_0_1) / M_PI * 180.0;

        // 计算起止两点切线的夹角
        if (!is_less_than(angle_0_1, at)) {
            // 创建一个新的AF_POINT类型几何点
            AF_POINT* nw = ACIS_NEW AF_POINT();
            // 设置该点的参数为当前起止段中点参数值
            nw->set_parameter(mid_t);
            // 设置新点的位置
            nw->set_position(eq.eval_position(mid_t));
            // 将新点插入当前点cur后
            nw->insert(cur, sense);
        }
        cur = nxt;
    }
}

void CheckSingularities(EDGE* e, AF_POINT* head, double dt, double at) {
    double sp = e->start_param(), ep = e->end_param();

    ENTITY_LIST faces;
    api_get_faces(e, faces);

    AF_POINT* cur = head;
    double nd = 0.0;

    bool flag = false;

    for (FACE* f = (FACE*)faces.first(); f != nullptr; f = (FACE*)faces.next()) {
        cur = head;

        if (f->geometry()->equation().type() == 1)
            continue;                                     // no singularity on plane
        else if (f->geometry()->equation().type() == 2) {  // cone
            SPAposition cone_singularity;
            cone* lc = (cone*)&f->geometry()->equation_for_update();
            cone_singularity = lc->get_apex();
            if (!lc->_IsCylinder) {
                SPAposition proj;
                SPAparameter guess = 0.0, actual = 0.0;
                e->geometry()->equation().point_perp(cone_singularity, proj, guess, actual);
                nd = (cone_singularity - proj).len();
                if (nd < 2.0 * dt) {
                    cur = head;
                    actual = e->geometry()->equation().param(proj);
                    while (cur->next(0) != head) {
                        if (is_less_than(cur->get_parameter(), actual) && is_greater_than(cur->next(0)->get_parameter(), actual)) {
                            AF_POINT* check_st = cur;
                            AF_POINT* check_ed = cur->next(0);
                            AF_POINT* nw = ACIS_NEW AF_POINT();
                            // 设置该点的参数为当前起止段中点参数值
                            nw->set_parameter(actual);
                            // 设置新点的位置
                            nw->set_position(cone_singularity);
                            // 将新点插入当前点cur后
                            nw->insert(cur, e->sense());

                            double check_dt = nd / 2.0;
                            break;
                        }
                        cur = cur->next(0);
                    }
                }
            }
        }
        else if (f->geometry()->equation().type() == sphere_type) {  // sphere
            sphere* localSphere = (sphere*)&f->geometry()->equation_for_update();
            SPAposition sphere_centre = localSphere->centre;
            SPAunit_vector sphere_pole_direction = localSphere->pole_dir;
            double sphere_radius = localSphere->radius;

            // singularities
            SPAposition upper_singularity = sphere_centre + sphere_radius * sphere_pole_direction;
            SPAposition lower_singularity = sphere_centre - sphere_radius * sphere_pole_direction;

            SPAposition proj;
            SPAparameter guess = 0.0, actual = 0.0;
            e->geometry()->equation().point_perp(upper_singularity, proj, guess, actual);
            nd = (upper_singularity - proj).len();
            if (nd < 2.0 * dt) {
                cur = head;
                actual = e->geometry()->equation().param(proj);

                while (cur->next(0) != head) {
                    if (cur->get_parameter() < actual && cur->next(0)->get_parameter() > actual) {
                        AF_POINT* check_st = cur;
                        AF_POINT* check_ed = cur->next(0);
                        AF_POINT* nw = ACIS_NEW AF_POINT();
                        // 设置该点的参数为当前起止段中点参数值
                        nw->set_parameter(actual);
                        // 设置新点的位置
                        nw->set_position(upper_singularity);
                        // 将新点插入当前点cur后
                        nw->insert(cur, e->sense());

                        double check_dt = nd / 2.0;
                        // CheckDistanceTolerance(e, check_dt, check_st, check_ed);
                        break;
                    }
                    cur = cur->next(0);
                }
            }
            cur = head;
            guess = 0.0, actual = 0.0;
            flag = false;
            e->geometry()->equation().point_perp(lower_singularity, proj, guess, actual);
            nd = (lower_singularity - proj).len();
            if (nd < 2.0 * dt) {
                cur = head;
                actual = e->geometry()->equation().param(proj);
                while (cur->next(0) != head) {
                    if (cur->get_parameter() < actual && cur->next(0)->get_parameter() > actual) {
                        AF_POINT* check_st = cur;
                        AF_POINT* check_ed = cur->next(0)->next(0);
                        AF_POINT* nw = ACIS_NEW AF_POINT();
                        // 设置该点的参数为当前起止段中点参数值
                        nw->set_parameter(actual);
                        // 设置新点的位置
                        nw->set_position(lower_singularity);
                        // 将新点插入当前点cur后
                        nw->insert(cur, e->sense());
                        // CheckDistanceTolerance(e, check_dt, check_st, check_ed);
                        break;
                    }
                    cur = cur->next(0);
                }
            }
        }
        else if (f->geometry()->equation().type() == torus_type) {  // torus
            SPAinterval torus_u_range = f->geometry()->equation().param_range_u();
            SPAinterval torus_v_range = f->geometry()->equation().param_range_v();
            double u_high = torus_u_range.end_pt(), u_low = torus_u_range.start_pt();
            double v_high = torus_v_range.end_pt(), v_low = torus_v_range.start_pt();
            SPAposition proj;
            SPAparameter guess = 0.0, actual = 0.0;
            flag = false;
            if (f->geometry()->equation().singular_u(u_high)) {  // singularities can only appear at end of param u range
                SPAposition upper_singularity = f->geometry()->equation().eval_position(SPApar_pos(u_high, (v_high + v_low) / 2.0));
                cur = head;
                e->geometry()->equation().point_perp(upper_singularity, proj, guess, actual);
                nd = (upper_singularity - proj).len();
                if (nd < 2.0 * dt) {
                    cur = head;
                    actual = e->geometry()->equation().param(proj);
                    while (cur->next(0) != head) {
                        if (cur->get_parameter() < actual && cur->next(0)->get_parameter() > actual) {
                            AF_POINT* check_st = cur;
                            AF_POINT* check_ed = cur->next(0)->next(0);
                            AF_POINT* nw = ACIS_NEW AF_POINT();
                            // 设置该点的参数为当前起止段中点参数值
                            nw->set_parameter(actual);
                            // 设置新点的位置
                            nw->set_position(upper_singularity);
                            // 将新点插入当前点cur后
                            nw->insert(cur, e->sense());
                            break;
                        }
                        cur = cur->next(0);
                    }
                }
            }
            if (f->geometry()->equation().singular_u(u_low)) {  // singularities can only appear at end of param u range
                SPAposition lower_singularity = f->geometry()->equation().eval_position(SPApar_pos(u_low, (v_high + v_low) / 2.0));
                cur = head;
                guess = 0.0, actual = 0.0;
                flag = false;
                e->geometry()->equation().point_perp(lower_singularity, proj, guess, actual);
                nd = (lower_singularity - proj).len();
                if (nd < 2.0 * dt) {
                    cur = head;
                    actual = e->geometry()->equation().param(proj);
                    while (cur->next(0) != head) {
                        if (cur->get_parameter() < actual && cur->next(0)->get_parameter() > actual) {
                            AF_POINT* check_st = cur;
                            AF_POINT* check_ed = cur->next(0)->next(0);
                            AF_POINT* nw = ACIS_NEW AF_POINT();
                            // 设置该点的参数为当前起止段中点参数值
                            nw->set_parameter(actual);
                            // 设置新点的位置
                            nw->set_position(lower_singularity);
                            // 将新点插入当前点cur后
                            nw->insert(cur, e->sense());
                            break;
                        }
                        cur = cur->next(0);
                    }
                }
            }
        }
        else {  // spline
            SPAinterval spline_u_range = f->geometry()->equation().param_range_u();
            SPAinterval spline_v_range = f->geometry()->equation().param_range_v();
            double u_high = spline_u_range.end_pt(), u_low = spline_u_range.start_pt();
            double v_high = spline_v_range.end_pt(), v_low = spline_v_range.start_pt();
            SPAposition proj;
            SPAparameter guess = 0.0, actual = 0.0;
            flag = false;
            if (f->geometry()->equation().singular_u(u_high)) {
                SPAposition singularity_1 = f->geometry()->equation().eval_position(SPApar_pos(u_high, (v_high + v_low) / 2.0));
                cur = head;
                guess = 0.0, actual = 0.0;
                e->geometry()->equation().point_perp(singularity_1, proj, guess, actual);
                nd = (singularity_1 - proj).len();
                if (nd < 2.0 * dt) {
                    cur = head;
                    actual = e->geometry()->equation().param(proj);
                    while (cur->next(0) != head) {
                        if (cur->get_parameter() < actual && cur->next(0)->get_parameter() > actual) {
                            AF_POINT* check_st = cur;
                            AF_POINT* check_ed = cur->next(0)->next(0);
                            AF_POINT* nw = ACIS_NEW AF_POINT();
                            // 设置该点的参数为当前起止段中点参数值
                            nw->set_parameter(actual);
                            // 设置新点的位置
                            nw->set_position(singularity_1);
                            // 将新点插入当前点cur后
                            nw->insert(cur, e->sense());
                            break;
                        }
                        cur = cur->next(0);
                    }
                }
            }
            if (f->geometry()->equation().singular_v(v_high)) {
                SPAposition singularity_2 = f->geometry()->equation().eval_position(SPApar_pos((u_high + u_low) / 2.0, v_high));
                cur = head;
                guess = 0.0, actual = 0.0;
                e->geometry()->equation().point_perp(singularity_2, proj, guess, actual);
                nd = (singularity_2 - proj).len();
                if (nd < 2.0 * dt) {
                    cur = head;
                    actual = e->geometry()->equation().param(proj);
                    while (cur->next(0) != head) {
                        if (cur->get_parameter() < actual && cur->next(0)->get_parameter() > actual) {
                            AF_POINT* check_st = cur;
                            AF_POINT* check_ed = cur->next(0)->next(0);
                            AF_POINT* nw = ACIS_NEW AF_POINT();
                            // 设置该点的参数为当前起止段中点参数值
                            nw->set_parameter(actual);
                            // 设置新点的位置
                            nw->set_position(singularity_2);
                            // 将新点插入当前点cur后
                            nw->insert(cur, e->sense());
                            break;
                        }
                        cur = cur->next(0);
                    }
                }
            }
            if (f->geometry()->equation().singular_u(u_low)) {
                SPAposition singularity_3 = f->geometry()->equation().eval_position(SPApar_pos(u_low, (v_high + v_low) / 2.0));
                cur = head;
                guess = 0.0, actual = 0.0;
                e->geometry()->equation().point_perp(singularity_3, proj, guess, actual);
                nd = (singularity_3 - proj).len();
                if (nd < 2.0 * dt) {
                    cur = head;
                    actual = e->geometry()->equation().param(proj);
                    while (cur->next(0) != head) {
                        if (cur->get_parameter() < actual && cur->next(0)->get_parameter() > actual) {
                            AF_POINT* check_st = cur;
                            AF_POINT* check_ed = cur->next(0)->next(0);
                            AF_POINT* nw = ACIS_NEW AF_POINT();
                            // 设置该点的参数为当前起止段中点参数值
                            nw->set_parameter(actual);
                            // 设置新点的位置
                            nw->set_position(singularity_3);
                            // 将新点插入当前点cur后
                            nw->insert(cur, e->sense());
                            break;
                        }
                        cur = cur->next(0);
                    }
                }
            }
            if (f->geometry()->equation().singular_v(v_low)) {
                SPAposition singularity_4 = f->geometry()->equation().eval_position(SPApar_pos((u_high + u_low) / 2.0, v_low));
                cur = head;
                guess = 0.0, actual = 0.0;
                e->geometry()->equation().point_perp(singularity_4, proj, guess, actual);
                nd = (singularity_4 - proj).len();
                if (nd < 2.0 * dt) {
                    cur = head;
                    actual = e->geometry()->equation().param(proj);
                    while (cur->next(0) != head) {
                        if (cur->get_parameter() < actual && cur->next(0)->get_parameter() > actual) {
                            AF_POINT* check_st = cur;
                            AF_POINT* check_ed = cur->next(0)->next(0);
                            AF_POINT* nw = ACIS_NEW AF_POINT();
                            // 设置该点的参数为当前起止段中点参数值
                            nw->set_parameter(actual);
                            // 设置新点的位置
                            nw->set_position(singularity_4);
                            // 将新点插入当前点cur后
                            nw->insert(cur, e->sense());
                            break;
                        }
                        cur = cur->next(0);
                    }
                }
            }
        }
    }
}

/*-----离散圆形曲线-----*/
/**
 * const ellipse* curve：指向ellipse对象的指针，代表一个圆形曲线（用的是椭圆类）
 * double st, double ed：椭圆参数形式的起始参数t和结束参数t
 * double dt：距离容差
 * double at：角度容差
 * int sense：生成点的方向 0为正向，1为反向
 * AF_POINT* head：指向第一个AF_POINT对象的指针，作为结果链表的头部
 */
void facet_edge_circle(const ellipse* curve, double st, double ed, double dt, double at, int sense, AF_POINT* head) {
    AF_POINT* cur = head;

    // 圆半径
    double radius = curve->GetMajorAxisLength();

    // 防止参数增量大于2Π
    double delta_t_max = std::fabs(ed - st) / 3.0;

    // 初始化当前t为起始参数t
    double cur_param = st;
    /*auto val = 1.0 - dt / radius;
    val = std::min(1.0, std::max(-1.0, val));
    auto angle = 2.0 * acos(val);
    angle = std::min(at * M_PI / 180.0, angle);*/
    // 获取圆弧的参数步进值
    double angle = GetSegAngle(dt, at, radius);
    // 如果是反向的就将角度取负，参数增量取绝对值较小的那个
    if (sense) {
        angle *= -1.0;
        angle = std::max(angle, -delta_t_max);
    }
    else {
        angle = std::min(angle, delta_t_max);
    }
    // 遍历曲线上的点直到超出结束点
    while ((cur_param < ed && sense == 0) || (cur_param > ed && sense == 1)) {
        if ((cur_param + angle >= ed && sense == 0) || (cur_param + angle <= ed && sense == 1)) {
            break;
        }
        cur_param += angle;

        // 创建一个新的AF_POINT对象
        AF_POINT* nw = ACIS_NEW AF_POINT();
        // 设置点的参数
        nw->set_parameter(cur_param);
        // 设置点的位置
        nw->set_position(curve->eval_position(cur_param));
        // 将该新点插入链表中
        nw->insert(cur, sense);
        // 当前节点的下一个节点变成当前节点
        cur = cur->next(0);
    }
}

/*-----离散椭圆曲线-----*/
/**
 * const ellipse* curve：指向ellipse对象的指针，代表一个椭圆曲线
 * double st, double ed：椭圆参数形式的起始参数t和结束参数t
 * double dt：距离容差
 * double at：角度容差
 * int sense：生成点的方向 0为正向，1为反向
 * AF_POINT* head：指向第一个AF_POINT对象的指针，作为结果链表的头部
 */
void facet_edge_ellipse(const ellipse* curve, double st, double ed, double dt, double at, int sense, AF_POINT* head) {
    AF_POINT* cur = head;

    // 获取椭圆线长短半轴长度
    double a = curve->GetMajorAxisLength();
    double b = a * curve->radius_ratio;

    // 防止参数增量大于2Π
    double delta_t_max = (ed - st) / 3.0;

    auto val_a = std::clamp(1.0 - dt / a, -1.0, 1.0);
    auto val_b = std::clamp(1.0 - dt / b, -1.0, 1.0);

    double delta_ta = 2.0 * acos(val_a);
    double delta_tb = 2.0 * acos(val_b);

    // 转换为弧度制
    double at_delta = at * M_PI / 180.0;

    auto compute_delta = [&](double t) -> double {
        double delta_t = 0;
        delta_t = 2 * acos(1 - dt / (a + (b - a) * abs(pow(sin(t), 3))));
        return delta_t;
        };

    // 从起始点开始
    double cur_param = st;
    while ((cur_param < ed && sense == 0) || (cur_param > ed && sense == 1)) {
        double t = cur_param;
        double delta_t = compute_delta(t);
        delta_t = std::min(delta_t, at_delta);
        // 考虑遍历方向，对步进值进行正负调整
        if (sense == 0)
            delta_t = std::min(delta_t, delta_t_max);
        else {
            delta_t *= -1.0;
            delta_t = std::max(delta_t, -delta_t_max);
        }
        // 从起始点遍历到终止点结束
        if ((cur_param + delta_t >= ed && sense == 0) || (cur_param + delta_t <= ed && sense == 1)) {
            break;
        }
        // 计算新点的参数位置
        cur_param += delta_t;

        // 创建一个新的AF_POINT对象
        AF_POINT* nw = ACIS_NEW AF_POINT();
        // 设置点的参数
        nw->set_parameter(cur_param);
        // 设置点的位置
        nw->set_position(curve->eval_position(cur_param));
        // 将该新点插入链表中
        nw->insert(cur, sense);
        // 当前节点的下一个节点变成当前节点
        cur = cur->next(0);
    }
}

/*-----离散螺旋曲线-----*/
/**
 * const helix* curve：指向helix对象的指针，代表一个螺旋曲线
 * double st, double ed：螺旋参数形式的起始参数t和结束参数t
 * double dt：距离容差
 * double at：角度容差
 * int sense：生成点的方向 0为正向，1为反向
 * AF_POINT* head：指向第一个AF_POINT对象的指针，作为结果链表的头部
 */

void facet_edge_helix(const helix* curve, double st, double ed, double dt, double at, int sense, AF_POINT* head) {
    AF_POINT* cur = head;

    // 获取helix的基本信息
    double taper = curve->taper();                // 螺旋曲线的锥度
    double par_scaling = curve->par_scaling();    // 螺旋曲线的参数缩放因子
    SPAvector start_disp = curve->start_disp();   // 距起始点的偏移量
    SPAunit_vector axis_dir = curve->axis_dir();  // 螺旋轴方向的单位向量

    // 防止参数增量大于2Π
    double delta_t_max = (ed - st) / 3.0;

    // 计算参数增量
    double cur_param = st;
    double r_len = (start_disp - (start_disp % axis_dir) * axis_dir).len();
    double rt = r_len - taper * (cur_param + M_PI) / (2.0 * M_PI * par_scaling);
    /*double val = 1.0 - dt / rt;
    val = std::min(1.0, std::max(-1.0, val));
    double delta_t = 2.0 * acos(val);
    delta_t = std::min(delta_t, at * M_PI / 180.0);*/
    double delta_t = GetSegAngle(dt, at, rt);
    // 考虑遍历方向，对步进值进行正负调整
    if (sense == 0)
        delta_t = std::min(delta_t, delta_t_max);
    else {
        delta_t *= -1.0;
        delta_t = std::max(delta_t, -delta_t_max);
    }
    while ((cur_param < ed && sense == 0) || (cur_param > ed && sense == 1)) {
        // 从起始点遍历到终止点结束
        if ((cur_param + delta_t >= ed && sense == 0) || (cur_param + delta_t <= ed && sense == 1)) {
            break;
        }
        // 计算新点的参数位置
        cur_param += delta_t;

        // 创建新点，并添加到链表
        AF_POINT* nw = ACIS_NEW AF_POINT();
        nw->set_parameter(cur_param);
        nw->set_position(curve->eval_position(cur_param));
        nw->insert(cur, sense);
        cur = cur->next(0);
    }
}

/*-----离散插值曲线-----*/
/**
 * EDGE* e：边缘对象
 * const helix* curve：指向intcurve对象的指针，代表一个插值曲线
 * double st, double ed：插值曲线参数形式的起始参数t和结束参数t
 * double dt：距离容差
 * double at：角度容差
 * int sense：生成点的方向 0为正向，1为反向
 * AF_POINT* head：指向第一个AF_POINT对象的指针，作为结果链表的头部
 */
void facet_edge_intcurve(const intcurve* curve, double st, double ed, double dt, double at, int sense, AF_POINT* head) {
    // 设置当前点为头节点
    AF_POINT* cur = head;
    // 获取插值曲线的三次B样条曲线
    bs3_curve bs3c = curve->cur();
    double delta_t = 0.0;
    // 获取B样条曲线的次数
    int curve_degree = bs3_curve_degree(bs3c);
    // 初始化节点数量
    int numkts = 0;
    // 初始化节点数组
    double* knots = nullptr;
    // 从曲线bs3c中获取节点信息，填充到numkts和knots中
    bs3_curve_knots(bs3c, numkts, knots);

    // 初始化参数为double最小值
    SPAparameter left_param = sense ? ed : st;

    SPAposition left_pos_3d = curve->eval_position(left_param);

    SPAparameter last_param = left_param;

    auto last_pos_3d = left_pos_3d;

    std::vector<double> wait_for_insert;

    int curve_begin_index, curve_end_index;

    if (sense) {
        curve_begin_index = std::upper_bound(knots, knots + numkts, ed) - knots;
        curve_end_index = std::lower_bound(knots, knots + numkts, st) - knots;
    }
    else {
        curve_begin_index = std::upper_bound(knots, knots + numkts, st) - knots;
        curve_end_index = std::lower_bound(knots, knots + numkts, ed) - knots;
    }

    for (int i = curve_begin_index; i < curve_end_index; i++) {
        SPAparameter right_param = knots[i];
        if (is_zero_mch(right_param - left_param))continue;
        auto right_pos_3d = curve->eval_position(right_param);
        // 节点去重
        SPAunit_vector vec_p0Tangent = curve->eval_direction(left_param, TRUE, TRUE);
        SPAunit_vector vec_p1Tangent = curve->eval_direction(right_param, TRUE, TRUE);
        // 计算起止两点切线的夹角
        double dot_0_1 = std::clamp(vec_p0Tangent % vec_p1Tangent, -1.0, 1.0);
        // 限制至[-1,1]
        double angle_0_1 = acos(dot_0_1) / M_PI * 180.0;
        auto dir_3d = normalise(right_pos_3d - left_pos_3d);
        if (ed < st) dir_3d = -dir_3d;
        auto a_1 = acos(std::clamp(vec_p0Tangent % dir_3d, -1.0, 1.0)) / M_PI * 180;
        auto a_2 = acos(std::clamp(vec_p1Tangent % dir_3d, -1.0, 1.0)) / M_PI * 180;
        double w = (a_2 + at / 2) / (a_1 + a_2 + at);
        // 获取两端点的参数域中点
        double mid_t = left_param + (right_param - left_param) * w;
        auto real_mid_3d = curve->eval_position(mid_t);
        auto mid_3d = left_pos_3d + (right_pos_3d - left_pos_3d) * w;

        if (dt * dt < (real_mid_3d - mid_3d).len_sq() || is_greater_than(angle_0_1, at)) {
            wait_for_insert.push_back(last_param);
            left_param = last_param;
            left_pos_3d = last_pos_3d;
        }
        last_param = right_param;
        last_pos_3d = right_pos_3d;
    }
    if (wait_for_insert.size() < curve_degree + 1) {
        int initial_split_segement_count = curve_degree + 1;
        delta_t = (ed - st) / (initial_split_segement_count);
        // if(sense) delta_t = -delta_t;
        for (int i = 1; i < initial_split_segement_count; i++) {
            double nw_t = st + i * delta_t;
            // 创建新点
            AF_POINT* nw = ACIS_NEW AF_POINT();
            // 设置新点的参数
            nw->set_parameter(nw_t);
            // 设置新点的位置
            nw->set_position(curve->eval_position(nw_t));
            // 将新点插入链表
            nw->insert(cur, sense);
            // 更新当前点到下一个点
            cur = cur->next(0);
        }
    }
    else {
        if (sense) std::reverse(wait_for_insert.begin(), wait_for_insert.end());
        for (auto x : wait_for_insert) {
            // 创建新点
            AF_POINT* nw = ACIS_NEW AF_POINT();
            // 设置新点的参数
            nw->set_parameter(x);
            // 设置新点的位置
            nw->set_position(curve->eval_position(x));
            // 将新点插入链表
            nw->insert(cur, sense);
            // 更新当前点到下一个点
            cur = cur->next(0);
        }
    }
    // 计算参数步进值
}